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Abstract 

This paper describes a novel numerical model aiming at solving moving-boundary 
problems such as free-surface flows or fluid-structure interaction. This model uses a 
moving-grid technique to solve the Navier-Stokes equations expressed in the arbi- 
trary Lagrangian-Eulerian kinematics. The discretization in space is based on the 
spectral element method. The coupling of the fluid equations and the moving-grid 
equations is essentially done through the conditions on the moving boundaries. Two- 
and three-dimensional simulations are presented: translation and rotation of a cylin- 
der in a fluid, and large-amplitude sloshing in a rectangular tank. The accuracy and 
robustness of the present numerical model is studied and discussed. 

Key words: Spectral element method, moving-boundary problem, ALE, 
moving-grid. 



1 Introduction 



With the advent of powerful computational resources like clusters of PCs 
or parallel computers the numericists are able to address more challenging 
problems involving multi-physics and multi-scale approaches. These problems 
cover a large spectrum of scientific and engineering applications. However, in 
this paper, for the sake of conciseness, we will restrict our attention to two 
specific problems, namely: free-surface flows and fluid-structure interaction. 
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Free-surface flows occur in many industrial applications: coating flows, vertical 
drawing of viscous fluids, jets, die flows, etc, and in environmental flows: ocean 
waves, off-shore engineering, coastal habitat and management, to name a few. 
Two review articles have been published in recent years and report the state- 
of-the-art of the field [1,2]. It can be observed that free-surface flows have 
been tackled by direct numerical simulation at low and moderate Reynolds 
numbers. This reality is essentially due to the nonlinear characters of the 
flow. On top of the nonlinearity associated to the Navier-Stokes equations 
themselves, here we deal with a complicated geometry which is changing in 
time and which is part of the solution itself. This accumulation of difficulties 
calls for elaborate algorithms and numerical techniques. 

Fluid-structure interaction has been recognized for a long time as a real chal- 
lenge. Indeed, this interaction is present in engineering problems like turbo- 
machinery, aerospace applications: buffeting, acoustics, and also in biomedical 
flows like blood flow in the coronary arteries. Fluid-structure interaction is also 
encountered in the field of vortex-induced vibrations having many important 
marine applications (e.g related to oil exploration, cable dynamics, etc.). It 
is only at the present time that this type of interaction for three-dimensional 
cases appears to be feasible as the necessary computing power becomes avail- 
able. On one hand, the computational fluid dynamics (CFD) codes integrate 
the full steady state or transient Navier-Stokes equations which govern the 
dynamics of a viscous Newtonian fluid. They mostly use finite volume or finite 
element approximations. On the other hand, the computational solid mechan- 
ics (GSM) codes integrate the dynamics of various solid models, incorporating 
for example, classical infinitesimal linear elasticity, nonlinear finite elasticity 
with large deformations, plasticity, visco-elasticity, etc. These problems are 
also highly nonlinear with respect to the complicated geometries at hand. The 
combination of the nonlinearities of the mathematical models for the con- 
stitutive relations and for the geometrical behaviour has called for a robust 
approach able to deal with all the complexities and intricacies. The finite el- 
ement method (FEM) with the isoparametric elements has emerged as the 
leading technology and methodology in GSM. 

In the present paper, the methodological framework is the same for the fluid 
and the solid parts and rests upon the spectral element method [3-6]. With this 
choice the space discretization is similar for both problems. As in free-surface 
flows and fluid-structure interaction the geometry is deforming and moving, it 
is needed to use the arbitrary Lagrangian-Eulerian (ALE) formulation [7-10]. 
This formulation allows to treat the full geometrical problem with respect to 
a reference configuration that is arbitrarily chosen. A mapping is introduced 
to ease the description of the current configuration with respect to a reference 
configuration. This process leads to an ALE velocity which will be related to 
a grid velocity. 
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In Section 2, the mathematical models will be presented with the associated 
weak formulations in the ALE context. Section 3 will be devoted to space and 
time discretizations. Section 4 will describe the numerical algorithms for the 
moving-grid technique. Section 5 will present numerical results and the final 
section will draw some conclusions. 



2 Mathematical model 

A moving boundary-fitted grid technique has been chosen to simulate the 
unsteady part of the boundary in our computations. In the particular cases 
dealt with in this paper, the unsteady part of the boundary can be either 
the free surface in case of free-surface flows, or for fluid-structure interaction 
problems, the interface between the fluid and the structure. This choice of a 
surface-tracking technique is primarily based on accuracy requirements. With 
this group of techniques, the grid is configured to conform to the shape of 
the interface, and thus adapts continuously — at each time step — to it and 
therefore provides an accurate description of the moving boundary to express 
the related kinematic and/or dynamic boundary conditions. 

The moving-boundary incompressible Newtonian fluid flows considered in this 
paper, are governed by the Navier-Stokes equations comprising the momen- 
tum equation and the divergence-free condition. In the ALE formulation, a 
mixed kinematic description is employed: a Lagrangian description of the mov- 
ing boundary, an Eulerian description of the fixed domain boundaries and a 
mixed description of the internal fluid domain. 

2. 1 The ALE kinematic framework 

This section will introduce the notations used in this paper to define the vari- 
ables and frames of reference related to the ALE formulation. The notations 
adopted hereafter are borrowed from [10, 11]. We denote by VLt the fiuid do- 
main subset of with = 2, 3 the space dimension, the subscript t referring 
to the time t as the fiuid domain is changing when its boundaries are mov- 
ing. Let us denote by VLq a reference configuration — for instance the domain 
configuration at initial time t = Iq. The system evolution is studied in the 
time interval / = [to,^]- The position of a point in the current fiuid domain 
Vtt is denoted by x — Eulerian coordinate — and in the reference frame VLq by 
Y — ALE coordinate. Let At be a family of mappings, which at each t & I 
associates a point Y G to a point x G fi^ 

At -.Vto^W^ ^VttC x(Y, t) = At{Y). (1) 
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At is assumed to be continuous and invertible on Qq and differentiable almost 
everywhere in /. The inverse of the mapping At is also continuous on Qq. 

The Jacobian matrix of the ALE mapping At is given by 

Ja = (2) 

and its determinant Jyi^ = det J_4j is the Jacobian of the mapping character- 
izing the metrics of Qt generated from the one of Qq. The Euler expansion 
formula gives the relationship between the Jacobian of the mapping At and 
the divergence of the ALE velocity w: 



dJAt 



dt 



J^Vx-w, \/{Y,t) eQoX I, (3) 



supplemented by the initial condition J^^ = 1 for t = t^. In real computations, 
w will be associated to the mesh velocity. The hyperbolic partial differential 
equation ([3]) highlights the important role played by the divergence of the 
mesh velocity in the time evolution of the mapping At- This particular point is 
emphasized in Section HI where one of the mesh-update techniques used in our 
simulations, enforces a divergence-free condition for w resulting in a constant 
in time Jacobian [12]. Furthermore Eq. ([3]) constitutes the evolution law 
for as in our formulation the ALE mesh velocity is calculated based on 
the essential boundary conditions of our problem in Qt, thereby defining the 
location of the grid nodes and the value of At- 

Considering a time-dependent scalar field / defined on fit x I, the notation 
df/dt\Y refers to the time derivative in the ALE frame, or in short the ALE 
time derivative expressed in Eulerian coordinates as opposed to the regular 
time derivative in Eulerian coordinates and denoted by df/dt\x- Finally, the 
ALE mesh velocity w is defined as 



w(x,t) 



at 



(4) 



It is worth noting that a standard application of the chain rule to the time 
derivative gives 



df 



Y dt 



+ w-Vx/. (5) 



dt 

The symbol Vx indicates the gradient operation applied to the scalar field / 
with respect to the Eulerian coordinate x. If w = 0, the mesh is fixed, and we 
recover the Eulerian description where d/dt\Y is the classical time derivative 
d/dt\y^. If w = u where u is the fluid velocity field, we obtain the Lagrangian 
description and d/dt\Y is the particle derivative. Eq. ([5]) allows to generalize 
the Reynolds transport theorem for a time-dependent volume integral of a 
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scalar field 

dt (Xt ) Int^dt ^ ) 



(6) 



2.2 The strong ALE formulation for the Navier-Stokes equations 



The governing equations of our moving-boundary problem in the ALE kine- 
matic description, for an incompressible Newtonian fluid flow occupying a fluid 
domain Qt whose boundary dQt is evolving with time, are the Navier-Stokes 
equations which in strong form and in the Eulerian kinematic description read 



du 

'dt 



u- Vxu = -Vx|)+2z/Vx-D,(u) + f, V(x,t) G X /, (7) 
Vx-u = 0, V(x,t)GfitX/, (8) 



where u is the velocity field, p the pressure field (normalized by the constant 
fluid density p and relative to zero ambient), u the kinematic viscosity of the 
fluid, Dx(u) = |(VxU + VxU-^) the rate-of-deformation tensor and f the body 
force per unit mass, with the superscript T indicating the transpose. Eq. ([7]) 
expresses the conservation of momentum and the divergence-free condition 
([H]) is the continuity equation in its simplified form for an incompressible flow. 
Equations ([7])-([H]) are valid in the internal fluid domain Qt, and have to be 
supplied with boundary conditions on the boundary dQt and the problem 
being unsteady, an initial condition is also required. The initial velocity field 
is chosen as 

u(x, t = to) = u°(x), Vx G r^j, with fit = fito = r^o, (9) 

such that Vx ■ u° = 0. Let the boundary dQt be split into two non-overlapping 
parts dQf = dQf U dQf. In the sequel we will consider the two following types 
of boundary conditions 

u = g(t), on dnf and Vt G /, (10) 
cr-n = -j9l-n + 2z/Dx(u)-n = h(t), on dnf and \/t e I , (11) 

where cr is the stress tensor, I the identity tensor and h the local unit outward 
normal vector to dflf Eq- (fTOl) is an essential boundary condition on dQf of 
Dirichlet type. In the cases of free-surface flows and fluid-structure interaction 
problems, which are of particular interest for us, dQf corresponds to fixed or 
prescribed moving solid walls where a no-slip condition has to be satisfied. 
Eq. (fTTj) is a natural boundary condition on dQf . For free-surface flows and 
fluid-structure interactions, dQf represents prescribed inflow and/or outflow 
depending on the situations considered, but primarily the free surface itself 
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or the interface between the fluid and the structure, where a mechanical equi- 
hbrium has to be enforced. Therefore (fTTj) is a dynamic boundary condition 
(DBC) expressing the continuity of the normal stress at the moving bound- 
ary. If free-surface flows are envisaged and if the surface tension is neglected, 
h = —poh where po is the pressure of the surrounding fluid, taken as zero in 
the sequel. 

Using Eq. ([5]), we can recast the strong form of the conservation of the mo- 
mentum of the Navier-Stokes equations in the ALE frame 



du 

'dt 



+ (u-w)-VxU= -Vxj5 + 2z/Vx-Dx(u) + f, V(x,t)GaxJ, (12) 



Y 



the divergence-free condition ([8]), the initial and boundary conditions (p!)- (fTT!) 
remaining unchanged in the ALE kinematic description. Indeed, boundary 
conditions are related to the problem and not to the kinematic description 
employed, be it Eulerian, Lagrangian or arbitrary Lagrangian-Eulerian. Nev- 
ertheless the ALE mesh velocity w is to a certain extent part of the unknown 
flelds of the problem even though some freedom in moving the mesh makes 
the ALE technique so attractive. The details related to the treatment and the 
computation of the mesh velocity are presented in Section |H 



2.3 The weak ALE formulation for the moving-boundary problem governed 
by the Navier-Stokes equations 



Based on the strong formulation of the moving-boundary problem described 
in Section 12. 2[ one can derive the more appropriate weak transient ALE for- 
mulation. In a standard approach, flrst are introduced the spaces of test and 
trial functions used to express the initial problem in its weak form based on 
the reference conflguration Qq- We may note that the spaces of test and trial 
functions for the pressure are identical as no essential Dirichlet condition is 
being imposed on this fleld. This space is the space of functions that are square 
Lebesgue-integrable on the domain Qt and is denoted by In general 

the velocity does not necessarily vanish on the domain boundary; in our par- 
ticular case the existence of a non-homogeneous essential Dirichlet boundary 
condition on dflf leads us to consider different spaces for the test and trial 
functions for the velocity fleld. The solution for the velocity u, of the problem 
fl7|)- f|TT|) will be searched for directly in the Sobolev space of trial functions 
H^{QtY deflned as follows 

H^^tf = {u e L^ntY, Vxu, GL2(a)^withz = l,...,ci, U|^^^=g}, 

(13) 

and corresponding to the current conflguration Qf The reference conflguration 
Qq will be used to build the velocity test functions v, which will be taken in 
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the space Hq-^^QoY with 

Hl^{no) = {veL\no), VJ^eL^QoY, ^Ux,=0}, (14) 



to satisfy a homogenous Dirichlet condition on dQ^. Subsequently, the dy- 
namics of the test functions on the configuration is obtained by means of 
the existing inverse of the mapping At- Therefore the velocity test functions 
on the configuration at time t are the set of functions (v o A^'^) with v be- 
longing to HQjy{floY. In the sequel, the notation (v o At'^) is kept in order 
to emphasize two key points. First, the treatment of the weak form of the 
time derivative du/dt\Y in Eq. (fT2l) is based on the essential property that v 
is not time dependent and consequently d-v/dt\Y = 0. Second, such formu- 
lation highlights the path to follow when practically implementing the weak 
ALE formulation. Indeed, the time dependency is fully incorporated in the 
inverse mapping A^^ and the functions v remains the same as the ones used 
in fixed-grid problems. Moreover, in a more general framework where At is 
still invertible but less regular, this formulation still holds and one only needs 
to care for the regularity of the functions v and not of the functions (vo^^^). 
With the notations and spaces introduced, the weak transient ALE formula- 
tion reads: 

Find (u(t),p(t)) G H^iVttY x L'^i^t) such that for almost every t > to 

/ (vo^-i) . ud(^+ / (v o ^-1) . Vx ■ [uu - uw] d(^ = 
at Jcit Jut 

f (pVx ■ (v o A;') - 2uB^{u) : Vx(v o A;')) dQ (15) 

+ [ f ■ (v o Ai^) dn+ I h ■ (v o At^) ddn, Vv e i7o^^(^o)^ 
Jfit Jon? 



and 

-/ gVx-udl] = 0, yqeL'^int). (16) 

The above set of equations has to be intended in the sense of distribution in 
the interval t > to, therefore justifying the qualifier "for almost every t > to" , 
see [13] for greater details. As expected the DBC ffTTl) appears 'naturally' in 
the weak formulation above, corresponding to the last term on the right-hand 
side of f[T^ and being the only 'surface term' as the spatial integration is lim- 
ited to dQ^. In addition, the DBC ([TT]) defines the reference pressure level and 
therefore no additional condition on the mean value of the pressure is required. 
Finally, it is well known (see [13] for instance) that a non-homogeneous Dirich- 
let boundary condition engenders a compatibility condition that the field u 
has to satisfy. The origin of this condition is that, in order to be compatible 
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with ([8]) the boundary condition (fTOl) imposes 

/ u-nddn= [ g{t)-hddn+ [ u{t) ■ nddil = 0, \/teI. (17) 
Jdnt Jdnf Jdnf 

Eq. (fT7|) is a consequence of (fT6|) with g = 1. 

In order to ease the discretization of the continuous weak equations (fT5|) - 
fll6l) . we introduce the following notations and bilinear forms, such as a scalar 
product defined by 

(u,v):=/ (vo^-i).udl], VvGi/oVW, (18) 

a so-called 'energy bilinear form' 

^(u,v) := 2u f D,(u) : Vx(vo^-i)dfi, Vv G //q't^W, (19) 

a bilinear form related to the weak incompressibility constraint 

i3(v,p) := - / pV, ■ (v o A;') dn, Vv G i/o'7?(^o)', (20) 

a trilinear form corresponding to the nonlinear convective term 

C(v; u, w) := / (v o A^^) ■ ■ [uu - uw] dO, Vv G Hl^iQoY, (21) 

and finally two linear forms, the first one being related to the source term f 

^(v) := / f ■ (v o A;') dn, Vv G Hl^inoY, (22) 

and the second one to the non-homogeneous natural boundary condition (fTTl) 
on dnf 

7^-(v) := / h- (vo^-i)dan, W^eH'iQoY. (23) 

With these notations, the continuous weak ALE form of our moving-boundary 
Navier-Stokes problem can be recast as 

Find (u(t),p(t)) G H\,{VltY x L'^{Vtt) such that for almost every t > to 

^(u,v)+^(u,v) + i3(v,p) + C(v;u,w)=^(v)+7^'^(v), Vv G ifoV(^o)", 

(24) 

i3(u,g) = 0, VgGL2(a). 

(25) 
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3 Numerical technique and discretization 

Moving-boundary problems, either free-surface or fluid-structure interaction, 
represent a real challenge not only for the mathematician but also for the 
numericists. As presented in Section |21 the weak formulation of the problem 
( I24l) -( l25i) is an evidence of its complexity. In this section, particular emphasis 
will be put on the numerical space discretization of this arduous continu- 
ous problem. The general case with non-homogeneous natural and essential 
Dirichlet boundary conditions is dealt with, justifying the authors' choice of 
a very detailed presentation. The particular case of steady problems with 
non-homogeneous Neumann conditions and homogeneous Dirichlet boundary 
conditions was first formulated by Ho and Patera in [14] in their study of 
free-surface flows dominated by inhomogeneous surface-tension effects. Fur- 
thermore, R0nquist extended the previous work of Ho and Patera to the more 
general case of steady free-surface flow problems with non-homogeneous Neu- 
mann and Dirichlet boundary conditions [15]. The specificities related to the 
treatment of transient problems is highlighted in the present paper, which to 
our knowledge is not available in the current literature. 



3.1 Spectral element discretization 

The first step in the spectral element method (SEM) discretization consists 
in subdividing the fluid domain fi^ = fi^ U dVtf into E non-overlapping ele- 
ments {^t}e=i- ■'■^ sequel we will assume that this elemental subdivision 
is maintained for all values of t in the interval J, therefore meaning that no 
re-meshing procedure is applied and leading to 



A re-meshing procedure for problem using SEM is presented in [12] and can 
be used if needed. Each element VL^ involves a mesh constructed as a ten- 
sor product of one-dimensional grids. Although each space direction may be 
discretized independently of the others, without loss of generality we will con- 
sider only meshes obtained with the same number of nodes in each direction, 
denoted by + 1, corresponding to the dimension of the space of A^th-order 
polynomials. To describe the discretization process accurately, we adopt the 
same procedure as in [5] and define the following spaces 




for e = 1, ■ ■ ■ , -E, 



(26) 



d 



Y := H^n^Y, 



(27) 
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3.2 Galerkin approximation 



We apply the Galerkin approximation to our Navier-Stokes problem in the 
ALE formulation in its weak form fl24p - (12S]) . and therefore select finite di- 
mensional polynomial subspaces X^r, Y/v and Zj^f to represent X, Y and Z 
respectively. A staggered- grid approach with element based on Fn — IPAr-2 
spaces for the velocity and pressure field respectively, allows to eliminate com- 
pletely the spurious pressure modes [16]. In this context, the finite dimensional 
functional spaces are defined as 

XN:=XnF%^E, YN:=Yn¥%^E, Z^ := Z nFN^2,E, (28) 

with 



Pm e = {010 € L'^i^t)', <P\n^ is a polynomial of degree < M, Ve = 1, ■ ■ ■ , E}, 

(29) 

where the superscript d in fl28p reflects the fact that test and trial velocity 
functions are d-dimensional. With these notations the Galerkin approximation 
of (IMD-dSSD reads 

Find (uAr(t),pAr(t)) eYnX Zm such that for almost every t >to 

^ (uat. Vat) + A{un, Vat) + B{vn, Pn) + 

C(vAr; Utv, Wat) = J^Ar(vAr) + 7^^(vAr), Vvat E X^, (30) 

B{uN,qN) = 0, VgAT e Zn, (31) 



with 



(uat. Vat) = V / uat ■ (vat o ^) dQ, Vva^ e Xat (32) 

e=l -^^t 
E 

-^iv(v7v) = E / ■ (v^ ° A') df^, Vv^ e Xjv (33) 

e=l ■'^t 

n'^.iyM) = E / h;v ■ (v;v o Ai^) ddn, Vv^ e Xn (34) 
Jan!'"' 



fAT and Wat being the projection of the fields f and w onto the finite dimen- 
sional space P^^;. 

The integrals within each of the spectral elements {^1)^=1 ^^"^ bound- 
aries are performed in a discrete manner using Gaussian quadra- 
ture rules. More specifically, all the terms in fl5Ul) - flHTl) are integrated using a 
Gauss-Lobatto-Legendre (GLL) quadrature rule except for the two terms in- 
volving the bilinear form B where functions discretized in PAr-2,£; appear. For 
these two terms, namely the pressure term and the divergence-free condition, 
a Gauss-Legendre (GL) quadrature rule is chosen. Therefore, the Pat — PAr-2 
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Navier-Stokes problem in the ALE form is finally stated as 
Find {uj\f{t),pi\f(t)) G Iat x Z^v such that for almost every t > to 

^ (uat, VAr)^ fj^^ + ^Ar,GLL(uAr, Vat) + BN,GLi'VN , Pn) + 

Cn,gll{^n] uat, Wat) = Tn,gll{^n) + 'H'^,gll{^n), Vvat G Xat, (35) 

BN,GL{viN,qN) = 0, VgAT G Zat. (36) 

To simplify the notations in the sequel, we will drop the subscript GLL and 
unless being explicitly specified, whenever an integration rule is required, the 
GLL one is implicitly being used. 



3.3 Semi-discrete Navier-Stokes moving-boundary problem in the ALE form 

In order to formulate the semi-discrete version of our moving-boundary prob- 
lem governed by the Navier-Stokes equations in the ALE form, we introduce 
the two tensor-product bases on the reference element Cl := [— 1, l]"^ and for 
the sake of simplicity we will choose the same discretization order in each 
space direction A^: 

• the Gauss-Lobatto-Legendre Lagrangian interpolation basis of degree A^ 

(37) 

to expand the velocity field discretized on the A^th-order GLL grid; 

• the Gauss-Legendre Lagrangian interpolation basis of degree N — 2 

(38) 

to expand the pressure field discretized on the GL grid of order N — 2. 

The expressions of the one-dimensional GLL and GL Lagrangian interpolant 
polynomials 7r(^) and ru{C) appearing in (P7I) - (I5S]) can be found in [5]. The 

polynomials {''Ti,j,k{^)}i^j^k=o ^^^1 {wij,A:(C)}ij7fc=i ^^^^ serve as bases for the 
functions in the spaces X^, Yjsf and Zn- 

N 

UN{x{^),t) = ^^jk{t) T^ijA^), V(^,t)G(^X/, (39) 

i,j,k=0 
N-1 

PN{^{C),t) = P^Ai) ^i,J,k{C), V(C,t)Gfix/, (40) 

i,j,k=l 
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where x = x*^ is the location of the point considered in the spectral element 
fif , {uijfc(t)}^- ;j^o the set of nodal values of the velocity field on the GLL grid 
of ilf and {Pijk{'t)}i^j~k=i the set of nodal values of the pressure field on the GL 
grid of Q^. It is important to note that the time-dependency of the discretized 
velocity uat and pressure pa? is not solely accounted by the time-dependent 
nodal values of these two fields. Indeed, due to the motion of the grid, the 
mapping between the position in the reference element Cl and the spectral 
element at time t is also time-dependent via the ALE mapping At- If we 
note Ail (resp. A^q) the mapping from the reference element Cl onto (resp. 
JIq), we can write 

x^ = X^(0, \/itt)eCixi, (41) 

Y^ = Ml{^), Witt) en X I, (42) 

x^ = A(Y^), V(Y^^)Gfigx/, (43) 

leading to following relationship between the different mappings 

= Ato Ml, Vt G /. (44) 

Eq. (jUj) shows that the second origin of the time-dependency of fl39|) and pUj) , 
after the one due to the set of GLL and GL nodal values, is the moving-grid 
technique via the time-dependency of the ALE mapping At- 

Before embarking on the final process of semi-descritizing the equations for 
the moving-boundary problem, a last issue needs to be addressed: the treat- 
ment of the non- homogeneous Dirichlet boundary condition (ITU]) on dQf . First 
of all and as mentioned earlier, the non- homogeneity of ffTOj) leads to differ- 
ent spaces for the trial and test functions for the velocity field, Xm and Yn 
respectively. Therefore the basis fl39|l developed for Xn is not suitable for 
the solution u^{t) of the problem fl5n]) - fl^ sought in Yjy. As presented ear- 
lier, the non- homogeneous Dirichlet boundary condition imposes to satisfy the 
compatibility condition (fTTl) whose discrete version reads 

E 

y unU) -hddn = 0, VtG/. (45) 

Let Ufo^Tv be a (piecewise) polynomial defined on the discrete boundary dVtl 
(e = 1, . . . , _E) and such that its nodal boundary values are equal to the bound- 
ary data g(t). In practice, the GLL Lagrangian interpolation bases defined on 
the element boundaries are used to expand VLb,N'i however, in the rest of the 
inner domain these functions are zero. By construction, Ub^jv satisfies (145|) . 
Setting mn = ^q,n + ^b,N, the problem reduces to finding uo,Ar in the space 
Yo^N ■■= Hljy{VttY n P^^^. Therefore the problem ([30D-([3l]) can be reformu- 
lated as follows 
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Find (uo,7v(t),PAr(t)) G Yq^n x Zn such that for almost every t > to 



— (uo,7V, Vtv)^ + ^Ar(uo,Ar, Vat) + Bn,Gl{^N, Pn) + 



dt 

Cn{^N] Uo,7V, Wtv) = J'7v(v7v) + H^ivN) + 1^i,n{^n, t) , Vvat G Xn, (46) 
BN,GL{'^0,N,qN) =1^2,N,GL{qN,t), V^Ar G Zat, (47) 



with 



d_ 

^iv(U6,Ar(t), Vat) - CAr(vAf; U6,Ar(^), Wat), Vvat G Xat, (48) 



T>l^Ni'VN,t) = -— (Ub,Af(t), VAr)Ar 



and 

1^2,N,GL{qN,'t) = -BN,GL{^b,N{t),qN), ^Qn^Zn. (49) 

The two time-dependent terms Vi^^ et V2^n,gl appearing in fH^ and (jUj) are 
due to the non-homogeneity of the Dirichlet boundary condition. These val- 
ues are related to the values of the discrete field Uf,^Ar(t), which as mentioned 
earlier, are obtained from the values of the field g(t) from (fTOj) . 



We can now expand the trial velocity uq^n and the trial pressure pn onto the 
GLL-GL bases like in (15^ and (HO]) respectively. In the sequel we will drop 
the subscript in uq^at, no confusion being possible as the non- homogeneous 
Dirichlet boundary conditions is already accounted for in (H6|) - p7|) . The semi- 
discrete equations derived from P6l) - p7|l are 

^(M(t)u(t)) = -K(t)u(t) - C(u(t), w(t),t)u(t) + D^(t)p(t) + Fi(t), (50) 
-D(t)u(t) = F2(t). (51) 



The matrices appearing in (l50l) - (l5Tl) are all time-dependent: M is the mass 
matrix, K the stiffness matrix, C the discrete convective operator involving 
the velocity field u and the ALE mesh velocity w, D-^ the discrete gradient 
operator and D the discrete divergence. and F_2 are two vectors accounting 
for the presence of the body force f and the time-dependent essential Dirichlet 
and natural non-homogeneous boundary conditions. 



3.4 Time discretization 



The set of semi-discrete equations ( !50|) - (!5T]) is discretized in time using finite- 
difference schemes in a decoupled approach. The computation of the linear 
Helmholtz problem — corresponding to the energy bilinear form A and the 
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stiffness matrix K — is integrated based on an implicit backward differentia- 
tion formula of order 2, the nonlinear convective term — corresponding to the 
trilinear form C and the matrix C — is integrated based on a relatively simple 
extrapolation method of order 2, introduced by Karniadakis et al. [17]. 

The moving-grid treatment requires the semi-discrete equations (|3U]) - flHT]) to 
be supplemented with an equation computing the mesh nodes update 

dx , , 

with X being a vector containing the ci- dimensional mesh nodes position at 
time t. The integration of Eq. fl52|) necessitates the knowledge of the values 
of the mesh velocity provided by the moving-grid technique employed. Two 
particular techniques are presented in detail in Section HI The time-integration 
of Eq. ( l52l) is based on an explicit and conditionally stable Adams-Bashforth 
of order 3. 



Lastly the treatment of the pressure relies on a generalized block LU decom- 
position with pressure correction [18,19]. 

The temporal order of the overall splitting scheme has proved to be equal to 
two for fixed-grid problems. The grid motion induces a limited reduction of 
the global temporal order, which is found to fall between 1.5 and 2 for the 
simulations presented in Section [5l The reasons for this reduction in the global 
order of the method is still not clearly understood. 



3. 5 Specificities pertaining to free-surface flows and fluid- structure interac- 
tion 



Up to this point, the treatment of our moving-boundary problem was kept 
to a level general enough to encompass both the free-surface flow and fluid- 
structure interaction problems. At this stage, it appears natural to provide 
the specificities pertaining to each of these two sub-problems. 

These specificities lie primarily in the natural and essential Dirichlet bound- 
ary conditions imposed to the system. For free-surface flows with no surface- 
tension effects — either of normal or tangential type — and with no inflow nor 
outflow — closed system, both the natural and essential Dirichlet boundary 
conditions become homogeneous — g = on dQf and h = on dQ^, lead- 
ing to a drastic simplification of the problem. More precisely, both vectors F^^ 
and F^2 vanish in the semi-discrete formulation (!50!) - (l5T!) of the problem. For 
fluid-structure interaction problems, the natural boundary condition on the 
interface between the flow and the structure is provided by the dynamics of the 
structure, that can be evaluated by the SEM and the Newmark method [20]. 
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In the sequel, we will focus our attention on the flow problem for both of these 
two cases of interest. 



4 Moving-grid techniques 

When considering moving-boundary problems tackled in an interface-tracking 
and ALE frame, the moving boundary dVt^ is treated in a Lagrangian way 
whereas the boundary OVt^ which is fixed or subject to a prescribed motion 
g(t), is studied in an Eulerian frame. Accordingly, such method allows large- 
amplitude motions of the moving boundary, by generating a shape-conformed 
grid. Hence, an accurate and simple application of the boundary conditions 
on dVtt is easily accessible. 

As a consequence of the ALE kinematics, the boundary conditions imposed 
on the mesh velocity w read 

w ■ n = u ■ fi, on ^VL^ and Vt G /, (53) 

w = g(t), on and Vt e /. (54) 

Eq. (l53l) is a kinematic boundary condition (KBC) on the moving boundary 
traducing that dVt^ is a material surface with no transfer of fluid across it in the 
Lagrangian perspective. Eq. (l54l) expresses a kinematic boundary condition of 
no-slip type on the boundary of the domain which is not free to move. Given 
fl53l) - fl54l) . it appears that the freedom left for the choice of w lies in the values 
of this field in the internal fluid domain VLt and also on the tangential values 
of w on the moving boundary dO!^ . 

The computation of the mesh velocity w in the internal fluid domain VLt is 
the corner-stone of the moving-grid technique developed in the framework of 
the ALE formulation. The values of the mesh velocity being prescribed on the 
boundary dVLt as expressed by equations flS^ - fl511) . the evaluation of w in VLt 
can be obtained as the solution of an elliptic equation: 

£xW = 0, V(x, t) G X /, (55) 

where £y^ represents any elliptic operator based on the Eulerian coordinates 
x. Such elliptic equation constitutes a classical choice for calculating the mesh 
velocity [21]. Two types of elliptic equations are envisaged in this paper. The 
first elliptic operator used is a Laplacian operator Ax, and will be used ex- 
tensively in the fluid-structure interaction simulations presented in Section O 
More details about the use of this specific operator for the computation of 
the mesh velocity can be found in [20]. The second strategy relies on the as- 
sumption that the motion of the mesh nodes is equivalent to a steady Stokes 
flow, corresponding physically to an incompressible and elastic motion of the 
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mesh. Therefore the problem for the mesh velocity becomes a boundary- value 
steady Stokes problem with the mesh velocity w satisfying a divergence-free 
condition 



Vx-w = V(x,t)GfitxJ. (56) 

The justifications of this additional constraint imposed to the mesh velocity 
problem is presented in detail in [12]. The free-surface flow simulations of a 
sloshing in a three-dimensional tank due to the gravity presented in Section 
15. 3[ were carried out using this second strategy for w. 



5 Numerical simulations and results 

In this section we will present numerical results for three problems: the steady 
Stokes equations in curved subdomains, the motion of a cylinder in a square 
cavity and the sloshing in a three-dimensional tank. 

5.1 Accuracy in curved domains 



We want to check the error evolution in the square domain Q = [—1,1]^ 
decomposed in curved subdomains (elements). To this aim, let us consider the 
steady Stokes equations 

- VxP + Axu + f = 0, V(x,t)Gax/, (57) 

Vx-u = 0, V(x,t)GfijxJ. (58) 

The exact solution is given by 

Ux = — cos {'JTx/2) sin {7iy /2) , (59) 

Uy = sin {nx/ 2) cos {ny/ 2) , (60) 

p = — vr sin (7rx/2) sin (7r?//2) , (61) 

when the body force term is chosen as 

fx = -ir^ cos (7rx/2) sin (7ry/2) , fy = 0. (62) 



Instead of using a regular square grid composed of elements with edges par- 
allel to the lines of the Cartesian axes, we performed the computation on the 
deformed mesh [22,23] displayed in Figured] (left). Contours of the norm of the 
velocity field for the computed solution of problem (l57j) - (!58|) , are presented on 
Figured] (right). Figure [2] shows the evolution of the relative error in if^-norm 
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for the velocity and in L^-norm for the pressure field, with respect to an in- 
creasing polynomial degree A^, for two cases — E = 2x2 elements and E = 4x4 
elements. The convergence is slower than the one obtained with a mesh divided 
in several regular square subdomains. First, we still achieve the exponential 
decrease of the relative error when the polynomial degree increases (which is 
typical of spectral or p-convergence [5]). Second, the convergence is faster when 
increasing the number of spectral elements E, as previously observed in [24] 
(which is equivalent to /i-convergence in finite-element terminology [5]). 
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Fig. 1. Square domain O with internally deformed subdomains with £' = 4x4 
spectral elements and = 20 (left) and the velocity magnitude (right) 
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Fig. 2. Relative error in i/^-norm for the velocity and in L^-norm for the pressure 
field 

The same computation has been carried out with a geometry Q' obtained by 
the transformation of coordinates of the unit square Q = [—1,1]^ with sine 
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functions (Fig. [3]) 

x' = X + a sin (vrx) sin (ny) , (63) 
y' = y + asm (ny) sin (vry) , (64) 

with {x,y) G ^2, {x',y') G and a = 1/10. Here, the deformation of the 
geometry not only involves the interior of the subdomains but also the domain 
boundaries. 

The remarks made for the first computation, corresponding to the square 
domain, are still relevant for this geometry. The same behavior of the con- 
vergence is obtained as one can observe on Figure |H The deformation of the 
boundaries induces obviously a slower convergence in comparison with the 
square domain. Nevertheless, the important result is that the spectral conver- 
gence is maintained (Fig. Hj) even with a deformation of the domain involving 
its boundaries, which is a mandatory feature when solving moving-boundary 
problems. 
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20 and = 30 (left) and the velocity magnitude 



5.2 Motion of a cylinder inside a square cavity 



We solve the Navier-Stokes equations ([7])-([HD in a two-dimensional square 
cavity. A schematic view of this cavity is given in Figure O with fixed exterior 
walls. A circular cylinder is immersed into the fluid and is moving with a 
prescribed velocity. Two types of prescribed motions are studied. In the first 
case, we consider a cylinder in horizontal translation from the center of the 
cavity with a constant velocity. Denoting the boundary of the cylinder as Fcyi, 



18 




^^2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 

polynomial degree 

Fig. 4. Relative error in H^-novm for the velocity and in L^-norm for the pressure 
field 

we prescribe 

w^lr^yi = w^|r,„ = 1, (65) 
%|rcyi = Wy\r,^^ = 0. (66) 

Denoting the exterior walls as Text, we have 

u|ro,t = w|r,,t = 0. (67) 

We solve a time-dependent problem in order to study the evolution of the 
fluid motion caused by the translation of the cylinder in the square domain 
Q = [—1, 1]^. The Reynolds number based on a unit reference length and a 
unit reference velocity is Re = l/u = 100. The diameter of the cylinder is 
D = 0.28. The time step At is fixed to 0.005. The discretization uses a total 
number of elements equal to E' = 64 and the polynomial degree is = 12 
in each of the two directions. We obtain an unsteady evolution of the fluid 
motion and we observe a deformation of the fluid mesh as pictured in Figure 
[6l where appears the flow configuration for t = 0.25, 0.5 and 0.7. If we keep on 
moving the cylinder closer to the right wall, the mesh deformation becomes 
too large. We have also focused our attention on the evaluation of an artificial 
"acceleration" defined as ||u„_|_i — UnW^'^/At, where u„ denotes the velocity 
field at the time level n. Figure [7| displays the previous expression and the 
L^-norm of the velocity. The "acceleration" does not vanish, which means the 
solution does not become steady-state. This can be expected since the cylinder 
is always in motion inside the cavity. Moreover the "acceleration" increases 
when the cylinder gets closer to the right wall. 
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Fig. 5. Geometry of the fluid domain with the immersed cylinder 
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Fig. 6. The velocity component Ux and the corresponding streamlines (black solid 
lines) around a moving cylinder, Re = 100, for t = 0.25 (left), t = 0.5 (center) and 
t = 0.7 (right) 

In the second case, the cylinder at the center of the cavity is subject to a 
constant counter-clockwise angular rotation uj = 1 such that 



(68) 
(69) 



In Figure [8], we have maintained E = 64 and changed the polynomial degree 
to = 10. We exhibit the flow configuration for t = 0.5, 1.25 and 2.0. Like in 
the previous example, we conclude the solution does not reach a steady state 
due to the motion of the cylinder (Fig. We have successfully tested these 
kinds of motion for large distortions of the fluid mesh. 



5.3 Sloshing in a three-dimensional tank 



To show the adaptability of our numerical model based on a moving-grid 
technique in the ALE frame, the analysis of large-amplitude sloshing in a three- 
dimensional rectangular tank has been carried out. The tank has a square-base 
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Fig. 7. ||u„_|_i — u„||j^2/At and ||u||^2 versus simulation time t 
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Fig. 8. The velocity component Ux and the corresponding streamlines (black solid 
lines) around a moving cylinder, Re = 100, for t = 0.50 (left), t = 1.25 (center) and 
t = 2.0 (right) 

section of dimensions L x L and a schematic view of the geometry of the 
problem is presented in Figure [TOl The free-surface position is measured from 
the bottom of the tank: z = L + h{x,y,t), where h{x,y,t) is the free-surface 
elevation measured from its equilibrium position z = L. The initial shape of 
the free surface is varying only with x and is given by the nonlinear theory for 
finite-amplitude standing waves 



h{x, t) = a{t) cos{kx) cos{ujt) 
ka^it) cos{2kx) 
2tanh(A;L) 



9, , 3 cos(2u;t) — tanh^(A;L) 1 
cos^ cut + ^ — — ) , 



(70) 



where the wave number is k = 27r/A, the wave length A = 2L, the initial 
wave amplitude a(t = 0) = L/5, and u = yJgktanh.{kL) corresponding to the 
dispersion relation. For the sake of simplicity, we have taken g = 27rA tanh(/i;L) 
which leads to a period T of the non-viscous and irrotational waves equal to 
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one. Finally, the Reynolds number is based on the reference velocity y/gL and 
is expressed as Re = L^J'gL/ v. 

Free surface at equilibrium 
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Fig. 10. Geometry of the simulation set-up and rectangular tank 

A no-slip condition is imposed to the velocity field at the bottom of the tank 
2; = and free-slip conditions on the side walls x = 0, a; = L and ?/ = 0, 
y = L, likewise for the mesh velocity field w which corresponds to Eq. (15^ 
in the context of this problem. The Dirichlet boundary condition on the free 
surface for the mesh velocity is given by w = u, which includes Eq. f l53|) and 
also an additional condition on the tangential values of w. The initial velocity 
field is the irrotational solution at the maximum displacement of a standing 
wave — corresponding to zero for all velocity components. When starting the 
simulation, the top surface is set free and allowed to evolve in response to 
the dynamic and kinematic boundary conditions (fTTl) and (153!) respectively. 
The nonlinearity of this problem is introduced by both boundary conditions, 
through the shape of the free surface in (l53i) and through the normal to the 



22 



free surface in ( ITT]) . The motion of the free surface physically corresponds to 
a transfer of energy between the potential energy — maximum at the initial 
time — and the kinetic energy, leading to a decaying oscillatory phenomenon. 

Computations are performed with = 3^ = 27 spectral elements and a 
polynomial degree iV = 9 in all three directions, leading to a mesh comprising 
28^ nodes. The time step is taken equal to 0.001 and the simulation duration 
is 25 time units — based on the unit period T — or 25'000 iterations for 7 values 
of the Reynolds number Re = 50, 100, 250, 500, 750, I'OOO and 1'500. 

Using an energetic argument. Lamb [25] derived the approximate damping of 
a free wave due to viscosity function of time 

a{t) = a(0)e-2'^'='*. (71) 

Figure [TT] displays the computed relative wave amplitude a(t)/a{0) with re- 
spect to the simulation time t/T. The excellent linear fits obtained for the 
seven values of the Reynolds number are in perfect agreement with the ex- 
ponential viscous damping. Equation fl71l) also shows that the relative wave 
amplitude is proportional to the kinematic viscosity, if plotted in y-\og scale. 
This second point is verified in Figure [T2| where a{t)/a{0) is plotted against 
the inverse of the Reynolds number which is by definition proportional to the 
kinematic viscosity u. 

Those results are evidences of the robustness and accuracy of our moving-grid 
technique in handling large-deformations for moving-boundary problems such 
as the one considered here. 
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Fig. 12. Relative amplitude a(t)/a(0) against the inverse of the Reynolds number 
6 Conclusions 



A numerical model for solving two- and three-dimensional moving-boundary 
problems such as free-surface flows or fluid-structure interaction is proposed. 
This model relies on a moving-grid technique to solve the Navier-Stokes 
equations expressed in the arbitrary Lagrangian-Eulerian kinematics and dis- 
cretized by the spectral element method. A detailed analysis of the continuous 
and discretized formulations of the general problem in the ALE frame, with 
non-homogeneous and unsteady boundary conditions is presented. Particular 
emphasis was put on the weak formulation and its semi-discrete counterpart. 
The moving-grid algorithm which is one of the key ingredient of our numer- 
ical model, is based on the computation of the ALE mesh velocity with the 
same accuracy and numerical technique as the fluid velocity. The coupling 
between the Navier-Stokes computation and the one for the mesh velocity 
is effective through the problem boundary conditions. It is noteworthy that 
the coupling in the interior Navier-Stokes computation is effective through the 
modified convective term which is induced by what is happening at the bound- 
aries. Three numerical test results are presented in the two particular cases of 
interest, namely fluid-structure interactions and free-surface flows. First the 
influence of the deformation of the grid on the accuracy of the numerical model 
is evaluated. In a second problem, two motions (translation and rotation) of 
a cylinder immersed in a fluid is computed. Lastly, large-amplitude sloshing 
in a three-dimensional tank is simulated. The results obtained are showing 
very good with the theoretical results when available, therefore leading to a 
validation of our numerical model. 
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